Explorando a base de dados da Covid-19¶
Começamos a explorar, aqui, a base de dados Coronavirus 2019-nCoV. Só um pontapé inicial. Vamos carregar as bases de dados e montar alguns gráficos. Neste processo, veremos a biblioteca pandas para o tratamento dos dados e as ferramentas gráficas matplotlib.pyplot, seaborn, e plotly.express.
Objetivos de aprendizado:
Ambiente kaggle.
Carregando dados no kaggle.
Trabalhando com pandas.
Gráficos com matplotlib.pyplot, pandas e plotly.
Crescimento exponencial.
Regressão linear e mínimos quadrados
Criando um kernel¶
Após criar uma conta no https://www.kaggle.com/, você pode criar um novo kernel de um notebook em python de duas maneiras:
Clicar no ícone
<>, na lateral esquerda, para acessar a página de notebooks e clicar no botão+ New Notebookno canto superior direito da página, para criar um novo notebook seu. Em seguida, selecionar o tipo de kernel. No caso, selecionar o tipo notebook e a linguagem python.Abrir um dos diversos kernels disponíveis no kaggle e clicar no botão
Copy and Editno canto superior direito da página, para criar uma cópia desse kernel na sua conta.
Você pode acessar os seus kernels clicando no ícone <> e na aba Your Work.
Gravando o kernel¶
Durante a edição do kernel, você pode fazer commits regulares, para salvar cada versão do kernel, à medida que você vai escrevendo mais.
Entre cada commit, você também tem a opção de deixar o kaggle salvar constantemente o rascunho. Isso pode, ou não, estar habilitado por default. Recomendo conferir e habilitar. Isso é feito abrindo a barra lateral direita (clicar na seta ao lado do botão Commit), se ela não estiver aberta, clicar em Settings e selecionar a opção On em Internet.
Infelizmente, não encontrei uma opção de "salvar" o kernel manualmente, mas, felizmente, com a opção internet habilitada, isso é feito com extrema frequência.
De qualquer forma, quando for dar uma pausa no trabalho, vale clicar na seta no canto superior esquerdo, para sair da edição (ou ameaçar sair), que o kaggle vai avisar se a última edição foi gravada ou não.
Funcionamento do kernel¶
Cada kernel funciona (quase) como um caderno Jupyter. Uma diferença é que não é permitindo importar dados diretamente da internet, e.g. o seguinte comando não funciona no kaggle, mas funciona localmente ou em outras plataformas de computação em nuvem.
pandas.read_csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv')
Veremos como resolver isso em seguida.
Origem e atualização dos dados¶
Os dados são carregados a partir de links montados no seu diretório de dados.
Esses links são criados de duas maneiras:
Ao criar um kernel a partir de outro kernel, os arquivos/links deste são automaticamente transferidos como links para o seu ambiente, no diretório
../input/ou/kernel/input/.Durante a edição do seu kernel, clicar em
Filena barra de menu superior e, em seguida, emAdd or upload data.
No caso deste kernel, os dados do covid-19 vêm do kernel Novel Corona Virus 2019 Dataset, que possui links para extrair os dados do repósitorio externo 2019 Novel Coronavirus COVID-19 (2019-nCoV) Data Repository by Johns Hopkins CSSE.
Os dados populacionais vêm do kernel Population Dataset Country-Wise, que possui um link para extrair esses dados da base disponibilizada em Worldometer: Countries in the world by population (2020).
Os links para dados externos são atualizados a cada commit do seu kernel.
Como dito acima, o kaggle não permite a extração direta de dados externos em células do kernel. Isso evita que os dados sejam constantemente carregados sempre que a célula é executada. Isso seria um desperdício de banda de tráfego de internet. Principalmente quando lembramos que muitos desses dados podem estar sendo usados por centenas de pessoas ao mesmo tempo (e.g. centenas de usuários trabalhando em uma mesma competição do kaggle, em cima dos mesmos dados).
Dessa forma, os dados externos são atualizados uma vez por dia pelo kaggle, armazenados em seus servidores e disponibilizados para os seus usuários diretamente como arquivos.
Isso reduz drasticamente o tráfego de dados na internet, acelera a execução de cada kernel e mantém os dados diariamente atualizados.
Importando as bibliotecas¶
Primeiro, importamos todas as bibliotecas usadas neste kernel.
# This Python 3 environment comes with many helpful analytics libraries installed
# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python
# For example, here's several helpful packages to load in
import os # system tools
import datetime as dt # date and time tools
import numpy as np # numerical package
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import matplotlib.pyplot as plt # graphics library
import seaborn as sns # fancy graphics library
import plotly.express as px # interactive graphics library
dt_string = dt.datetime.now().strftime("%d/%b/%Y")
print(f"Atualização mais recente do kernel: {dt_string}")
Carregando os dados¶
Arquivos de dados do Covid-19¶
Criamos este kernel a partir do kernel Novel Corona Virus 2019 Dataset, clicando no botão New Notebook. Dessa forma, este kernel foi automaticamente criado com os dados daquele kernel acessíveis no diretório /kaggle/input/ (ou ../input).
Arquivo de dados populacionais¶
Também adicionamos a base de dados do kernel Population Dataset Country-Wise. Para adicionar essa base de dados, é necessário clicar em File na barra de menu, selecionar Add or upload data, digitar "Population Dataset Country-Wise" no campo de buscas e em seguida clicar nesta base quando ela aparecer na lista de resultados da busca. Após isso, um link para essa base de dados é automaticamente criado no diretório /kaggle/input/ deste kernel.
Visualizando o diretório dos arquivos disponíveis¶
Podemos ver todos esses arquivos do diretório /kaggle/input com o script a seguir.
# Input data files are available in the "../input/" directory (or "/kaggle/input/").
# Running this (by clicking run or pressing Shift+Enter) will list all files under the input directory:
#for dirname, _, filenames in os.walk('/kaggle/input'):
# for filename in filenames:
# print(os.path.join(dirname, filename))
Carregando os arquivos¶
Carregamos os dados como DataFrames do pandas, usando a função pandas.read_csv().
Utilizamos, também, o método rename(), dos DataFrames, para renomear as colunas dessas tabelas. A opção inplace=True faz com que o próprio DataFrame seja modificado, caso contrário ele não é modificado, apenas retorna um outro DataFrame modificado que pode ser associado a uma nova variável.
def get_url_jh(filename):
"""
Retorna o endereço URL do arquivo 'filename' no repositório da John Hopkins University no github
"""
return 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/' + filename
df = pd.read_csv(get_url_jh('time_series_covid19_confirmed_global.csv'))
#df = pd.read_csv('/kaggle/input/novel-corona-virus-2019-dataset/covid_19_data.csv',parse_dates=['Last Update'])
df.rename(columns={'ObservationDate':'Date', 'Country/Region':'Country'}, inplace=True)
df_confirmed = pd.read_csv(get_url_jh('time_series_covid19_confirmed_global.csv'))
df_recovered = pd.read_csv(get_url_jh('time_series_covid19_recovered_global.csv'))
df_deaths = pd.read_csv(get_url_jh('time_series_covid19_deaths_global.csv'))
df_confirmed.rename(columns={'Country/Region':'Country'}, inplace=True)
df_recovered.rename(columns={'Country/Region':'Country'}, inplace=True)
df_deaths.rename(columns={'Country/Region':'Country'}, inplace=True)
O banco de dados da população do países inclui dados de vários anos e, por algum motivo, estão duplicados. Pegamos apenas o ano 2020 e descartamos as duplicatas.
df_pop = pd.read_csv("../input/population-dataset-countrywise/population_dataset.csv")
df_pop_2020 = df_pop[df_pop['year'] == '2020'].drop('year', axis=1).drop_duplicates(subset='country_name', keep='first')
df_pop_2020.rename(columns={'country_name':'Country', 'population':'Population'}, inplace=True)
Visualizando o conteúdo dos arquivos¶
Visualizamos o início ou o final de cada arquivo de dados da seguinte forma:
df_confirmed.head()
df_recovered.head()
df_deaths.tail()
df_pop_2020
Objetos do pandas¶
Cada objeto acima é do tipo DataFrame, cada coluna ou linha é do tipo Series e os elementos podem ser de diversos tipos.
print(type(df_confirmed))
print(type(df_confirmed['Country']))
print(type(df_confirmed.iloc[0]))
print(type(df_confirmed.iloc[0,0]))
print(type(df_confirmed.iloc[0,1]))
print(type(df_confirmed.iloc[0,2]))
print(type(df_confirmed.iloc[0,4]))
Dados do Brasil¶
Podemos extrair os dados do Brasil selecionando-os pela coluna do país.
Vamos aproveitar e juntar os dados (confirmados, recuperados, decessos) do Brasil em um mesmo DataFrame, com uma coluna para cada tipo de caso. Descartamos, ainda, as primeiras colunas, para ficar apenas com os dados diários da epidemia.
Selecionando a linha com os dados do Brasil¶
Em pandas, podemos selecionar uma (ou mais) linha(s) de acordo com alguma condição booleana. No caso, queremos que a coluna 'Country' seja igual a 'Brazil'. Isso é feito com a condição
df_confirmed['Country'] == 'Brazil'
que irá retornar um objeto pandas.core.series.Series com os valores True ou False em cada linha, dependendo do resultado do teste.
print(type(df_confirmed['Country'] == 'Brazil'))
print(df_confirmed['Country'] == 'Brazil')
Utilizando essa condição na seleção das linhas, que é feita com o método loc do DataFrame, nos dá a linha desejada, que, no caso, é a de número 28:
df_confirmed.loc[df_confirmed['Country'] == 'Brazil']
Observe que as primeiras colunas contém informações do país. Apenas a partir da quarta linha é que temos os dados da epidemia. Podemos selecionar isso pelo método iloc[0,4:]. O primeiro argumento 0 indica que é a linha 0 (e única) do DataFrame com a linha do Brasil.
O resultado é um objeto do tipo Series, já que selecionamos uma única linha.
print(type(df_confirmed.loc[df_confirmed['Country'] == 'Brazil'].iloc[0,4:]))
print(df_confirmed.loc[df_confirmed['Country'] == 'Brazil'].iloc[0,4:])
Juntando os dados do Brasil¶
Podemos fazer isso para cada tipo de caso e montar um DataFrame com todos os dados da epidemia no Brazil:
df_Brazil = pd.DataFrame({'Confirmed':df_confirmed.loc[df_confirmed['Country'] == 'Brazil'].iloc[0,4:],
'Recovered':df_recovered.loc[df_confirmed['Country'] == 'Brazil'].iloc[0,4:],
'Deaths':df_deaths.loc[df_confirmed['Country'] == 'Brazil'].iloc[0,4:]
})
df_Brazil.tail()
Dados relativos¶
Se quisermos ver os dados relativos à população total, em porcentagem, podemos dividir os dados acima pela população total e multiplicar por 100.
Para acessarmos a população do Brasil em 2020 no DataFrame df_pop_2020, selecionamos a linha 'Brazil' da coluna 'Country':
df_pop_2020[df_pop_2020['Country'] == 'Brazil']
Para obtermos o valor da população, selecionamos 'Brazil' na coluna 'Population':
df_pop_2020[df_pop_2020['Country'] == 'Brazil']['Population']
Mas este é um objeto do tipo Series:
print(type(df_pop_2020[df_pop_2020['Country'] == 'Brazil']['Population']))
print(df_pop_2020[df_pop_2020['Country'] == 'Brazil']['Population'])
Não podemos simplesmente dividir os dados do Covid-19 por um objeto desse tipo. Para isso, podemos ler o seu conteúdo a partir da sua posição no DataFrame:
df_pop_2020[df_pop_2020['Country'] == 'Brazil'].iloc[0,1]
Mas isso é uma string.
Não podemos convertê-la diretamente para um inteiro por conta das vírgulas.
Primeiro devemos eliminar as vírgulas e em seguida transformar essa nova string em um inteiro:
int(df_pop_2020[df_pop_2020['Country'] == 'Brazil'].iloc[0,1].replace(',',''))
Agora podemos obter as porcentagens de casos 'confirmados', 'recuperado' e 'decessos'.
df_Brazil_rel = 100000*df_Brazil/int(df_pop_2020[df_pop_2020['Country'] == 'Brazil'].iloc[0,1].replace(',',''))
df_Brazil_rel
Selecionando parte dos dados de acordo com um limiar¶
Como a epidemia que se iniciou na China demorou a chegar no Brasil, assim como em outros países, é interessante, visando análises estatísticas, restringir os dados ao início da epidemia. Podemos fazer isso através do método ge(n), que retorna um novo DataFrame, com True ou False, dependendo do valor original ser maior ou igual do que o argumento n. Usamos isso como condição para obtermos os dados desejados.
Como os números de casos confirmados, recuperados e de decessos não devem ser decrescentes, não deveríamos precisar nos preocupar em selecionar o primeiro índice em que o valor deles é maior ou igual a n e usar esse índice para selecionar todos os dados a partir desse índice. Basta selecionarmos todos as linhas em que esse valor é maior ou igual a n.
Assim, para selecionarmos apenas do dados de casos confirmados no Brasil a partir do momento em que o número desses casos chegou a dez, fazemos
df_Brazil[df_Brazil['Confirmed'].gt(10)]
Mas pode haver erros nos dados, então o ideal seria descobrir o primeiro dia em que a condição ocorre e selecionarmos todos os índices desse em diante. Mas não vou me preocupar com isso neste momento.
Assim, temos
df_Brazil[df_Brazil['Confirmed'].gt(10)]
Função para extrair os dados de cada país¶
Ao invés de repetirmos o procedimento acima para cada país que desejemos analisar, podemos facilitar a extração desses dados definindo uma função apropriada.
Mas devemos tomar cuidado. Alguns países aparecem diversas vezes no DataFrame, separados por Province/State. Por exemplo, vejamos o Canadá:
df_confirmed.loc[df_confirmed['Country'].str.contains("Canada")]
Nesse caso, para obtermos o resultado total, devemos somar as colunas. Neste processo, descartamos as quatro colunas iniciais, que não são de dados numéricos das ocorrências.
df_confirmed.loc[df_confirmed['Country'].str.contains("Canada")].sum().iloc[4:]
Agora, sim, estamos aptos a definir a função.
Incluímos, nessa função, uma opção para selecionar o país; uma opção para retornar os dados absolutos (número de casos) ou relativos (porcentagens em relação à população total, ou em relação a outros números, digamos, a cada mil habitantes, etc.); e uma opção para selecionar os dados a partir do momento em que um determinado tipo de caso ultrapassar for maior ou igual a um determinado valor.
def get_country_data(country, percount = False, case = '', threshold = 0):
"""
Retorna um *DataFrame* do pandas com os dados conjuntos de um mesmo país.
Junta os dados dos três *DataFrames* de casos confirmados, recuperados e decessos
em um mesmo *DataFrame*.
Opção para retornar dados absolutos ou relativos à população total, em relação a
um número de `percount`, por exemplo `percount=`100` para porcentagem, ou
`percount=1000` para cada mil habitantes.
Opção para descartar os dados em que a quantidades de casos do tipo `case` é
menor do que o valor `threshold`.
"""
data = pd.DataFrame({'Confirmed':df_confirmed.loc[df_confirmed['Country'].str.contains(country)].sum(0).iloc[4:],
'Recovered':df_recovered.loc[df_recovered['Country'].str.contains(country)].sum(0).iloc[4:],
'Deaths':df_deaths.loc[df_deaths['Country'].str.contains(country)].sum(0).iloc[4:]
})
if threshold and case in ('Confirmed', 'Recovered', 'Deaths'):
data = data[data[case].ge(threshold)]
if percount:
data = percount*data/int(df_pop_2020[df_pop_2020['Country']==country].iloc[0,1].replace(',',''))
return data
Extraindo os dados da Itália¶
Com essa função, podemos facilmente extrair os últimos dados de países quaisquer.
Por exemplo, a seguir, os dados da Itália em números absolutos:
data = get_country_data('Italy')
data
Em relação a cada mil habitantes:
get_country_data('Italy', 100000).tail()
E a partir do momento em que o número de casos confirmados chega a 1000:
get_country_data('Italy', case='Confirmed', threshold=1000)
for country in ('Italy', 'Brazil'):
get_country_data(country).plot(figsize=(12,4), title=country)
plt.show()
Embelezando o gráfico¶
Só para embelezar um pouco o gráfico, podemos usar o estilo de gráfico do pacote seaborn. Escolhemos o estilo darkgrid:
sns.set_style("darkgrid")
A partir disso, montamos o gráfico como antes.
for country in ('Italy', 'US'):
get_country_data(country).plot(figsize=(12,4), title=country)
plt.show()
Gráficos interativos¶
Podemos construir gráficos interativos com diversos pacotes. Aqui usamos o pacote plotly, mais precisamente o módulo plotly.express. Com ele, podemos ampliar, reduzir ou selecionar uma parte do gráfico, visualiar os valores exatas passando o cursor em cima das linhas, etc.
Para traçar múltiplas linhas em um mesmo gráfico usando o plotly, é preciso alterar o DataFrame para ter apenas três colunas, uma para a abscissa, outra para identificar a linha e outra para a ordenada. No nosso caso, as linhas são tipos de casos ("confirmados", "recuperados" e "decessos"). Isso é obtido com o método pandas.melt().
Vejamos o seu efeito no DataFrame dos dados do Brazil.
wide_df = get_country_data('Brazil')
tidy_df = pd.melt(wide_df.reset_index(), id_vars='index')
tidy_df.rename(columns={'index':'Data', 'variable': 'Situação', 'value': 'Número de casos'},
inplace=True)
tidy_df
for country in ('Italy', 'US'):
wide_df = get_country_data(country)
tidy_df = pd.melt(wide_df.reset_index(), id_vars='index')
tidy_df.rename(columns={'index':'Data', 'variable': 'Situação', 'value': 'Número de casos'},
inplace=True)
fig = px.line(tidy_df, x="Data", y="Número de casos", color="Situação",
title = {'text': country, 'x':0.5, 'y':0.85, 'xanchor': 'center'},
width=800, height=400)
fig.show()
Podemos, ainda, fazer uma animação, com os parâmetros adequados.
# https://plotly.com/python-api-reference/generated/plotly.express.line.html
country = 'Italy'
wide_df = get_country_data(country)
tidy_df = pd.melt(wide_df.reset_index(), id_vars='index')
tidy_df.rename(columns={'index':'Data', 'variable': 'Situação', 'value': 'Número de casos'},
inplace=True)
px.line(tidy_df, x="Data", y="Número de casos", color="Situação",
title = {'text': country, 'x':0.5, 'y':0.85, 'xanchor': 'center'},
hover_name='Situação',
animation_frame="Data", animation_group='Situação',
range_x=['1/22/20','3/29/20'], range_y=[-10000,100000],
width=800, height=500)
Usando subgráficos¶
Subgráficos podem ser feitos de diversas maneiras. Abaixo montamos uma delas, usando o função subplots() do matplotlib.pyplot.
Observe que cada plot é feito como um método de um DataFrames, que funciona porque os gráficos do pandas são feitos em cima do matplotlib.pyplot.
fig, axes = plt.subplots(figsize=(12,12), nrows=3, ncols=2, sharey='row')
for num, country in enumerate(['Italy', 'US', 'France', 'Spain', 'Argentina', 'Brazil']):
get_country_data(country).plot(ax=axes[num // 2, num % 2], title=country)
plt.show()
Múltiplos gráficos em uma mesma figura¶
Podemos combinar os dados de vários países em uma mesma figura. Uma maneira é criar um DataFrame com todos os dados que queremos combinar. Uma outra maneira mais simples, combinando os dados de poucos países e aproveitando o que já fizemos, é a seguinte.
fig, ax = plt.subplots(figsize=(12,4))
get_country_data('Brazil')['Confirmed'].plot(ax=ax)
get_country_data('Italy')['Confirmed'].plot(ax=ax)
get_country_data('US')['Confirmed'].plot(ax=ax)
get_country_data('France')['Confirmed'].plot(ax=ax)
get_country_data('Spain')['Confirmed'].plot(ax=ax)
get_country_data('Argentina')['Confirmed'].plot(ax=ax)
plt.legend(['Brasil', 'Itália', 'EUA', 'França','Espanha', 'Argentina'])
plt.title('Casos confirmados')
plt.show()
Se quisermos combinar os vários tipos de casos, podemos fazer como a seguir.
fig, (ax1, ax2, ax3) = plt.subplots(figsize=(12,15), nrows=3, ncols=1)
get_country_data('Brazil')['Confirmed'].plot(ax=ax1, title='Confirmados')
get_country_data('Italy')['Confirmed'].plot(ax=ax1)
get_country_data('US')['Confirmed'].plot(ax=ax1)
get_country_data('France')['Confirmed'].plot(ax=ax1)
get_country_data('Spain')['Confirmed'].plot(ax=ax1)
get_country_data('Argentina')['Confirmed'].plot(ax=ax1)
ax1.legend(['Brasil', 'Itália', 'EUA', 'França','Espanha', 'Argentina'])
get_country_data('Brazil')['Recovered'].plot(ax=ax2, title='Recuperados')
get_country_data('Italy')['Recovered'].plot(ax=ax2)
get_country_data('US')['Recovered'].plot(ax=ax2)
get_country_data('France')['Recovered'].plot(ax=ax2)
get_country_data('Spain')['Recovered'].plot(ax=ax2)
get_country_data('Argentina')['Recovered'].plot(ax=ax2)
ax2.legend(['Brasil', 'Itália', 'EUA', 'França','Espanha', 'Argentina'])
get_country_data('Brazil')['Deaths'].plot(ax=ax3, title='Decessos')
get_country_data('Italy')['Deaths'].plot(ax=ax3)
get_country_data('US')['Deaths'].plot(ax=ax3)
get_country_data('France')['Deaths'].plot(ax=ax3)
get_country_data('Spain')['Deaths'].plot(ax=ax3)
get_country_data('Argentina')['Deaths'].plot(ax=ax3)
ax3.legend(['Brasil', 'Itália', 'EUA', 'França','Espanha', 'Argentina'])
plt.show()
Crescimento exponencial¶
No início de uma epidemia, o crescimento do número de casos tem uma forma exponencial.
Digamos, por exemplo, que
$$ q(t) = be^{at}, $$onde $q$ é uma determinada quantidade (e.g. número de casos confirmados, ou recuperados ou decessos), $t$ é a variável temporal em uma determinada unidade de tempo (digamos, de um dia) e $a$ e $b$ são parâmetros relacionados à epidemia.
Ao determinarmos os parâmetros $a$ e $b$, podemos estimar o valor da quantidade $q$ em um tempo futuro e estimar como essa quantidade varia a cada unidade de tempo.
Observe que o parâmetro $a$ tem dimensão física inversa do tempo e que $b$ tem a mesma dimensão física da quantidade $q$. A dimensão de uma quantidade costuma ser denotada por um colchete. Assim,
$$ [a] = \frac{1}{\mathrm{tempo}} = \frac{1}{T}, \qquad [b] = [q]. $$A exponencial, o logaritmo e as funções trigonométricas são funções adimensionais (no domínio e no contradomínio).
Como $a$ tem dimensão inversa do tempo e $t$ tem dimensão de tempo, vemos que o produto $at$ que aparece na exponencial é adimensional, fazendo sentido.
O que o parâmetro $a$ nos diz¶
A partir do parâmetro $a$, assumindo-o constante em um determinado intervalo de tempo, e sendo a solução da forma
$$ q(t) = be^{at}, $$obtemos as sequintes informações:
- A cada intervalo de tempo $\Delta t$, a quantidade aumenta para
ou seja, ela aumenta pelo fator
$$ e^{a\Delta t} $$- Dito de outra forma, a quantidade aumenta em
- O tempo de duplicação é dado por $t_2$ tal que
ou seja,
$$ t_2 = \frac{1}{a}\ln 2. $$- Se $e^{am\Delta t} = n$, com $n, m\in \mathbb{N}$, então
então a quantidade é multiplicada por $n$ a cada intervalo de tempo $m\Delta t$.
- Se $t_2$ duplica, então, obviamente, $3t_2/2$ triplica, $2t_2$ quadriplica, e assim por diante.
Taxa de variação e taxa de variação específica¶
A taxa de variação dessa quantidade em relação ao tempo é
$$ \frac{\text{d} q}{\text{d} t} = a q $$Dividindo isso por $q$, temos a taxa de variação específica, que é exatamente o fator exponencial $a$:
$$ \frac{\frac{\text{d} q}{\text{d} t}}{q} = a $$Obs: Em termodinâmica, quantidades assim são chamadas intensivas, pois não dependem da quantidade total do "material", e.g. densidade, temperatura. Já $q' = aq$, que depende da quantidade total $q$, é o que em termodinâmica é chamado de quantidade extensiva, e.g. massa, volume e massa térmica.
Aproximações¶
Ao determinarmos o parâmetro $a$, vemos que, a cada intervalo de tempo $\Delta t$, a quantidade $q$ apresenta uma variação da ordem de
$$ \Delta q \approx a q \Delta t $$e a quantidade atualizada fica da ordem de
$$ q + \Delta q \approx q + aq\Delta t = (1 + a \Delta t) q. $$Escrevendo $q_n = q(t_0+n\Delta t)$, $n\in \mathbb{N}$, vemos que
$$ q_{n+1} = (1+a\Delta t) q_n $$Temos, também
$$ \frac{q(t + 2\Delta t) - q(t + \Delta t)}{q(t + \Delta t) - q(t)} \approx \frac{a q(t+\Delta t) \Delta t}{a q(t) \Delta t} = \frac{q(t+\Delta t)}{q(t)} \approx \frac{(1 + a\Delta t) q(t)}{q(t)} = (1 + a\Delta t). $$Fazendo $(\Delta q)_n = q_{n+1} - q_n$, podemos escrever, mais geralmente,
$$ (\Delta q)_{n+1} \approx (1 + a\Delta t) (\Delta q)_n. $$De outra forma, temos
$$ \frac{(\Delta q)_{n+1}}{(\Delta q)_n} \approx (1 + a\Delta t). $$Ou seja, durante um crescimento exponencial, a variação a cada unidade de tempo também cresce exponencialmente:
$$ (\Delta q)_n \approx (1 + a\Delta t)^n (\Delta q)_0 $$Quanto maior $a$, mais dramático o crescimento.
Obs: Fixando $t>0$ e escolhendo uma sequência $\delta_n = t/n$, vemos que
$$ (1+a\delta_n)^n = \left(1 + \frac{a t}{n}\right)^n \rightarrow e^{at}, \quad n\rightarrow \infty. $$Crescimento exponencial variável¶
Na prática, a taxa de crescimento exponencial varia com o tempo,
$$ \frac{\text{d} q}{\text{d} t} = a(t) q $$o que nos leva a
$$ q(t) = b e^{\int_0^t a(s) \text{d} s}. $$Fases¶
Em relação à terminologia, uma endemia é quando o número de casos é essencialmente constante (ou pelo menos periódico, no caso de ser sazonal). Um exemplo é a gripe comum, conhecida mais tecnicamente como influenza.
Chamamos de epidemia quando o número de casos de uma determinada doença é crescente. Esse crescimento é tipicamente exponencial.
No início de uma epidemia, a taxa de crescimento $a(t)$ do número de casos começa pequena, mas é positiva e crescente. Nesse caso,
$$ (\Delta q)_{n+1} > (\Delta q)_n, $$e o número de casos diários (novos casos) também aumenta.
A taxa $a(t)$ continua crescendo até atingir um valor máximo, representando o crescimento mais rápido possível. Após esse ponto de inflexão, o número de casos continua crescendo, mas a uma taxa menor. Nessa fase,
$$ (\Delta q)_{n+1} < (\Delta q)_n, $$e o número de casos diários (novos casos) diminui.
O pico da epidemia ocorre quando $a(t) = 0$. Após esse ponto, $a(t) < 0$ e o número de casos totais também começa a diminuir,
$$ \frac{\text{d} q}{\text{d} t} = a(t) q < 0. $$Deste ponto em diante, dizemos que chegamos ao fim da epidemia.
A tendência é que ela se estabilize para uma endemia, mas, em certos casos, pode haver novas ondas de epidemia.
Aproximações¶
Localmente, ainda podemos escrever
$q(t+\Delta t) \approx (1 + a(t)\Delta t) q(t)$.
$\Delta q \approx a(t)q(t)\Delta t$.
Tempo de duplicação a partir de um instante $t$ é aproximadamente $\ln 2/a(t)$.
Escrevendo $a_n = a(t+n\Delta t)$, temos
$q_n \approx (1+a_n\Delta t)q_n$.
$(\Delta q)_n \approx a_nq_n\Delta t$.
$(t_2)_n \approx \ln/a_n$.
Crescimento exponencial e gráficos em escala logarítmica¶
- Quando uma função apresenta um crescimento exponencial, como $q=be^{at}$, é útil fazer um gráfico em escala logarítmica no eixo das ordenadas. Isso significa fazermos uma mudança de variáveis
de forma que, $ Q = \ln q = \ln (be^{at}) = \ln b + \ln( e^{at}) = \ln b + at,$ ou seja,
$$ Q = B + at, $$onde $B=\ln b$. Ou seja, o gráfico se torna uma reta e a inclinação da reta nos diz o fator $a$ de crescimento exponencial.
- Outra forma de verificar isso é plotando a taxa de variação $q'$ em função da quantidade $q$, já que
Nesse caso, temos retas que passam pela origem e não precisamos de um parâmetro análogo a $B$.
- Além disso, usando escala logarítmica nos eixos de $q'$ e $q$, temos
Nesse caso, a inclinação é sempre $1$ e o parâmetro $a$ entra apenas para transladar a reta verticalmente.
Obs: Escalas logarítmicas nos dois eixos também são úteis em leis de potência $x = t^p$, com $p$ entrando como sendo a inclinação da reta obtida.
Crescimento da epidemia do covid-19¶
Vamos visualizar o crescimento dos dados da epidemia do covid-19 e ver se conseguimos verificar se o crescimento é exponencial.
Para isso, vamos usar a escala logarítmica na ordenada, do número de casos confirmados, em relação ao tempo.
for country in ('France', 'Italy', 'US', 'Brazil'):
get_country_data(country, False, 'Confirmed', 10).plot(figsize=(12,5),
logy=True,
title=country)
plt.show()
Ou em versão interativa, com o plotly.express. É um tanto irritante que os argumentos do plotly sejam um pouco diferentes do matplotlib.pyplot (e.g. logy vs log_y, title, etc.), mas paciência.
for country in ('France', 'Italy', 'US', 'Brazil'):
wide_df = get_country_data(country, False, 'Confirmed', 10)
tidy_df = pd.melt(wide_df.reset_index(), id_vars='index')
tidy_df.rename(columns={'index':'Data', 'variable': 'Situação', 'value': 'Número de casos'},
inplace=True)
fig = px.line(tidy_df, x="Data", y="Número de casos", color="Situação", log_y=True,
title = {'text': country, 'x':0.5, 'y':0.85, 'xanchor': 'center'},
width=800, height=400)
fig.show()
Crescimento no Brasil¶
No momento, vamos usar o "olhômetro" para encontrar parâmetros adequados. Posteriormente, vamos usar regressão linear para encontrar os parâmetros $a$ e $b$ que melhor se ajustam aos dados, em um determinado intervalo.
Para efeito de ilustração, usamos os dados de casos "Confirmados", e a partir do momento em que esse número chega a dez. Selecionamos esses dados através da função get_country_data():
df_Brazil_confirmed_10 = get_country_data('Brazil', False, 'Confirmed', 10)['Confirmed']
df_Brazil_confirmed_10
Este é um objeto pandas do tipo Series:
type(df_Brazil_confirmed_10)
df_Brazil_confirmed_10.index
Escolhemos $b$ como sendo o primeiro valor dessa série,
$$ b = \mathtt{df}\_\mathtt{Brazil}\_\mathtt{confirmed}\_\mathtt{10[0]} = 13 $$e escolhemos $a$ (por tentativa e erro) como sendo
$$ a = 0.27 $$Essa escolha de $a$ nos dá
$$ e^a \approx 1.29693 $$Com ele, montamos a nossa aproximação exponencial
$$ be^{at} $$Fazemos isso construindo uma lista em python com esses valores diários e transformando essa lista em um objeto Series do pandas, usando o mesmo índice dos dados do Brasil.
Os índices são acessados pelo atributo index:
df_Brazil_confirmed_10.index
Assim, montamos a aproximação exponencial.
a = 0.27
b = df_Brazil_confirmed_10[0]
approx_data = [b*np.exp(a*j) for j in range(df_Brazil_confirmed_10.count())]
df_approx = pd.Series(approx_data, index = df_Brazil_confirmed_10.index)
df_approx
Agora, podemos montar o gráfico. Dessa vez, para ficar mais claro que os nosso dados são discretos (diários), plotamos os dados reais com pontos.
fig, ax = plt.subplots(figsize=(12,5))
df_Brazil_confirmed_10.plot(ax=ax, style='o', logy=True)
df_approx.plot(ax=ax, style='--')
plt.title('Casos confirmados no Brasil')
ax.legend(['casos confirmados', f'be^(at), b = {b:.1f}, a = {a:.2f}'])
plt.show()
Vemos, pelo gráfico, que $a=0.27$ é uma boa aproximação. Calculando
np.exp(a) - 1
vemos que o número de casos diários aumenta em aproximadamente 31%.
Calculando
np.log(2)/a
vemos que o número de casos diários duplica aproximadamente a cada dois dias e meio.
Dito de outra forma, ele quadriplica a cada cinco dias.
Mas observe, no gráfico, que, na última semana, esse crescimento tem sido menor. Um ajuste apenas nesses últimos dados iria nos revelar um fator $a$ menor.
Mas não devemos nos ater a esses dados como sendo muito precisos. Há vários problemas práticos de falta de kits de testes, casos claros de infectados que não são testados, e atraso no resultado dos testes.
Regressão linear¶
Vamos, agora, usar regressão linear para encontrar os parâmetros $a$ e $b$ que melhor se ajustam aos dados. Vamos trabalhar em cima dos dados de números confirmados no Brasil.
Regressão linear trata da busca de um modelo linear que melhor se ajusta a valores dados a~priori. Isso alcançado através do problema de mínimos quadrados.
Mínimos quadrados¶
Considere um problema linear da forma
$$ A\mathbf{u} = \mathbf{c}. $$Suponha que este não tem solução, mas que estamos interessados em achar a "melhor solução possível" $\mathbf{u}$, no sentido de diminuir a norma Euclidiana (ou distância) entre $A\mathbf{u}$ e $\mathbf{c}$, i.e.
$$ \min_{\mathbf{u}\in \mathbb{R}^2} \|A\mathbf{u} - \mathbf{c}\|_2, $$Isso é equivalente a minimizar o quadrado da distância, que é uma função contínua, coerciva e continuamente diferenciável, dada por
$$ f(\mathbf{u}) = \|A\mathbf{u} - \mathbf{c}\|_2^2 = \|A\mathbf{u}\|_2^2 - 2\langle A\mathbf{u},\mathbf{c}\rangle + \|\mathbf{c}\|_2^2. $$Os seus pontos de mínimo podem ser obtidos procurando-se os pontos críticos da função.
O gradiente dessa função é dado por
$$ \mathbf{\nabla} f(\mathbf{u}) = 2A^tA\mathbf{u} - A^t \mathbf{c}, $$onde $A^t$ é a transposta da matriz $A$.
Assim, os seus pontos criticos são dados pelas soluções do sistema
$$ 2A^tA\mathbf{u} = A^t \mathbf{c}. $$Tipicamente, $A$ é uma matriz com várias linhas (correspondendo aos vários dados disponíveis) e poucas colunas (correspondendo aos poucos parâmetros que queremos usar para ajustar o modelo). Por isso é comum não termos solução para o problema inicial.
Por sua vez, $A^tA$ é uma matriz quadrada e esperamos encontrar solução, tipicamente uma única solução.
Este sistema pode até ser degenerado e nos dar uma infinidade de mínimos, mas pelo menos uma solução iremos encontrar. Além disso, a função $f$ não tem pontos de máximo, a menos que $A$ seja nula.
Aproximação linear¶
No caso de buscarmos uma aproximação da forma
$$ Q = B + at $$no plano $t \times Q$, com $Q = \ln q$, para os dados $(t_j,Q_j)$ no logaritmo $Q_j=\ln q_j$ do número $q_j$ de casos confirmados a cada dia $t_j$ no Brasil, queremos resolver aproximadamente o sistema
$$ \begin{cases} B + at_1 & = & Q_1, \\ B + at_2 & = & Q_2, \\ \;\;\;\vdots & \vdots & \vdots \\ B + at_n & = & Q_n \end{cases} $$Em forma de sistema, isso significa resolver o problema linear com
$$ A = \left[ \begin{matrix} t_1 & 1 \\ \vdots & 1 \\ t_n & 1 \end{matrix}\right], \qquad \mathbf{u} = \left( \begin{matrix} B \\ a \end{matrix}\right), \qquad \mathbf{c} = \left( \begin{matrix} Q_1 \\ \vdots \\ Q_n \end{matrix} \right). $$A matriz $A$ é uma matriz de Vandermonte, montada a partir dos instantes em que as medições são fornecidas. As incógnitas são os parâmetros de ajuste. E o lado direito são os (logaritmos) dos números de casos.
Resolvendo numericamente¶
Para resolver o problema de mínimos quadrados numericamente, vamos usar a função np.linalg.lstsq() do numpy.
Para isso, devemos montar a matriz $A$ e o vetor $\mathbf{c}$ acima.
Primeiro transformamos os dados de df_Brazil_confirmed_10, que é um objeto do tipo Series, em um array do numpy, usando o método to_numpy(), e calculamos o logaritmo dos seus coeficientes, para obter $\mathbf{c}$
c = np.log(df_Brazil_confirmed_10.to_numpy(dtype=int))
c
Em seguida, criamos o vetor de tempos. Como os dados são diários, i.e. $t_j = j-1$, simplesmente.
t = np.array([j for j in range(len(c))])
t
Agora, montamos a matriz de Vandermonde.
A = np.vstack([t,np.ones(len(t))]).T
A
Finalmente, resolvemos via
a_r, B_r = np.linalg.lstsq(A, c, rcond=None)[0]
print(a_r, B_r)
b_r = np.exp(B_r)
print(a_r, b_r)
Observe que isso ficou bem próximo do obtido "a olho".
Transformando em uma objeto do pandas, temos
regression_data = [b_r*np.exp(a_r*j) for j in range(df_Brazil_confirmed_10.count())]
df_regression = pd.Series(regression_data, index = df_Brazil_confirmed_10.index)
df_regression
df_regression - df_approx
Fazendo o gráfico, vemos
fig, ax = plt.subplots(figsize=(12,5))
df_Brazil_confirmed_10.plot(ax=ax, style='o', logy=True)
df_approx.plot(ax=ax, style='--')
df_regression.plot(ax=ax, style=':')
plt.title('Casos confirmados no Brasil')
ax.legend(['casos confirmados', 'aproximação via olhômetro', 'regressão linear'])
plt.show()
Erro quadrático e resíduos¶
O erro quadrático, ou soma dos quadrados dos resíduos, de um problema de mínimos quadrados é o valor de $\|A\mathbf{u} - \mathbf{c}\|_2^2$ para a melhor aproximação $\mathbf{u}$ encontrada, ou seja, é a soma do quadrado dos resíduos, que são os erros de cada medição,
$$ r_j = (A\mathbf{x})_j - y_j. $$Vamos comparar os erros quadráticos de cada aproximação:
SS_res_approx = np.linalg.norm(np.log(df_approx.to_numpy(dtype=int))
- np.log(df_Brazil_confirmed_10.to_numpy(dtype=int)))**2
SS_res_approx
SS_regression_res = np.linalg.norm(np.log(df_regression.to_numpy(dtype=int))
- np.log(df_Brazil_confirmed_10.to_numpy(dtype=int)))**2
SS_regression_res
Obs: No caso da regressão, temos a opção de obter a soma dos quadrados dos resíduos SS_res diretamente da função np.linalg.lstsq()
Erro quadrático médio¶
Como o erro quadrático é a soma dos resíduos, ele tende a aumentar a medida em que adicionamos dados na amostras/medição, apesar da aproximação ficar potencialmente melhor. Por conta disso, é costume se considerar o erro quadrático médio, obtido dividindo-se o erro quadrático pelo número de amostras.
Assim, dadas aproximações $\tilde Q_j$ para os valores $Q_j$ da amostra, independente delas terem sido obtidas via mínimos quadrados ou outra modelagem, temos que
- o resíduo de cada dado é a diferença
entre o valor $\tilde Q_i$ do modelo e o valor $Q_i$ do dado, para cada abscissa $t_i$, das $N$ amostras $i=1, \ldots, N$.
- A soma dos quadrado dos resíduos é o total
- O erro quadrático médio toma a média dos $r_i^2$:
R quadrado¶
O erro quadrático e o erro quadrático médio, no entanto, são valores absolutos e, por si só, não são boas medidas de quão bom o modelo é.
Para isso, uma primeira medida da qualidade do método é o R quadrado, ou coeficiente de determinação, definido por
$$ R^2 = 1 - \frac{\mathrm{SS}_\mathrm{res}}{\mathrm{SS}_\mathrm{tot}}, $$onde $\mathrm{SS}_\mathrm{tot}$ é a variância dos dados,
$$ \mathrm{SS}_\mathrm{tot} = \sum_i (Q_i - \bar Q)^2, $$e $\bar Q$ é a média dos valores:
$$ \bar Q = \frac{1}{N}\sum_i Q_i. $$Tipicamente, $R^2$ é um valor entre $0$ e $1$. E quanto mais perto de $1$, melhor a aproximação.
Uma interpretação para $R^2$ é a de indicar a fração dos dados que é bem "explicada", ou "corretamente prevista", pelo modelo.
Vamos calcular esse valor no nosso caso.
Q_medio = np.mean(np.log(df_Brazil_confirmed_10.to_numpy(dtype=int)))
SS_tot = np.linalg.norm(np.log(df_Brazil_confirmed_10.to_numpy(dtype=int)) - Q_medio)**2
R_sq = 1 - SS_regression_res/SS_tot
R_sq
Por fim, é costume traçar o gráfico da modelagem exibindo o valor de R quadrado:
fig, ax = plt.subplots(figsize=(12,5))
df_Brazil_confirmed_10.plot(ax=ax, style='o', logy=True)
df_regression.plot(ax=ax, style=':')
plt.title('Casos confirmados no Brasil')
ax.legend(['casos confirmados', f'regressão linear ($R^2 = {R_sq:.3f}$)'])
plt.show()
Em coordenadas lineares¶
Observe, no entanto, que em coordenadas lineares (reais), a aproximação exponencial não é tão boa:
fig, ax = plt.subplots(figsize=(12,5))
df_Brazil_confirmed_10.plot(ax=ax, style='o')
df_regression.plot(ax=ax, style=':')
plt.title('Casos confirmados no Brasil')
ax.legend(['casos confirmados', f'regressão linear ($R^2 = {R_sq:.3f}$)'])
plt.show()
Isso pode se dar por uma série de motivos.
Primeiro que, na realidade, o parâmetro de crescimento específico $a=a(t)$ não é constante. O crescimento está mais próximo de um crescimento logístico.
Segundo que o número de casos confirmados não corresponde ao número de casos de infectados.
Terceiro que essa diferença tende a aumentar a medida em que o número de casos aumenta e isso pressiona o sistema de saúde a reduzir, proporcionalmente, o número de testes e a aumentar o tempo até se obter o resultado. (O meu, por exemplo, demorou 15 dias, quando normalmente deveria durar entre dois e três dias).
Quarto que, em escala logarítmica, o error para valores grandes da quantidade é reduzido, diminuindo, relativamente, o peso deles no ajuste.
Exercícios¶
Faça um gráfico do número de casos da China, usando a função
get_country_data(). Determine o que há de errado no gráfico e implemente uma solução.Modifique a função
get_country_data, em relação aos argumentoscaseethreshold, para selecionar os dados a partir do primeiro índice em que o valor mínimo do caso selecionado é alcançado, para evitar que erros nos dados façam algumas datas futuras serem descartadas, causando um "buraco" na série de dados.Modifique a função
get_country_data()para incluir um novo argumento denominadolocalque selecione a Província/Estado (Province/State) do país. Faça esse argumento default para a string nula "" e a busca pela Província/Estado só ter efeito caso este argumento não seja "". Exemplifique isso com algumas províncias do Canadá.Modifique a função
get_country_dataincluindo um argumento denominadotranscorridospara que as datas sejam substituídas por dias transcorridos desde a data inicial, portanto começando com dia 0 e continuando com 1, 2, 3, etc. Faça esse argumento terFalsecomo valor default e que a substituição ocorra apenas se ela forTrue.Usando a modificação acima, faça um gráfico conjunto dos casos confirmados de vários países (a sua escolha, pelo menos dez), onde a abscissa é o número de dias transcorridos desde que o número de casos ficou maior ou igual a dez.
Repita o gráfico anterior, mas em escala logarítmica no eixo das ordenadas.
Sejam $\Delta t = 1\,\textrm{dia}$, $t_0$ o primeiro dia em que o número de casos confirmados é pelo menos 10, $q_n = q(t_0 + n\Delta t)$ o número de casos confirmados a cada $n$-ésimo dia após $t_0$, e $\Delta q_n = q_{n+1} - q_n$, $n\in \mathbb{N}$. Faça o gráfico de $\Delta q_n$ por $q_n$ (ou por $q_{n+1}$, se preferir), tanto em escala linear (em ambas as variáveis) como em escala logarítmica (em ambas as variáveis), para pelo menos dez países a sua escolha. De preferência que tenham características diferentes. Faça comentários.
Considere a equação $q'(t) = a(t) q(t)$ no caso em que $a=a(t)$ é variável.
a. Ache uma aproximação para a relação $(\Delta q)_{n+1}/(\Delta q)_n$ envolvendo apenas $a_n$, $a_{n+1}$ e $\Delta t$, onde $a_n = a(n\Delta t)$ e $a_n'= a'(n\Delta t)$.
b. Em seguida, usando a equação $q' = a q$, escreva a segunda derivada $q''$ em função de $a'$, $a$ e $q$.
c. Compare a condição $(\Delta q)_{n+1}/(\Delta q)_n < 1$ com a condição $q''< 0$. Comente.
Observações¶
Vale, em um segundo momento, verificar a base Coronavirus 2019-nCoV, que "conserta" (linhas diferentes relacionadas ao mesmo país, com nomes escritos de maneira diferentes) e completa (latitude e longitude) alguns dados de um dos arquivos da base acima, a saber csse_covid_19_daily_reports.csv.